Monitoring, detection, and surveillance system using principal component analysis with machine and sensor data

ABSTRACT

A computing device updates an estimate of one or more principal components for a next observation vector. An initial observation matrix is defined with first observation vectors. A number of the first observation vectors is a predefined window length. Each observation vector of the first observation vectors includes a plurality of values. A principal components decomposition is computed using the initial observation matrix. The principal components decomposition includes a sparse noise vector s, a first singular value decomposition vector U, and a second singular value decomposition vector ν for each observation vector of the first observation vectors. A rank r is determined based on the principal components decomposition. A next principal components decomposition is computed for a next observation vector using the determined rank r. The next principal components decomposition is output for the next observation vector and monitored to determine a status of a physical object.

CROSS-REFERENCE TO RELATED APPLICATIONS

The present application claims the benefit of 35 U.S.C. § 119(e) to U.S.Provisional Patent Application No. 62/462,291 filed on Feb. 22, 2017,the entire contents of which is hereby incorporated by reference.

BACKGROUND

Low-rank subspace models are powerful tools to analyze high-dimensionaldata from dynamical systems. Applications include tracking in radar andsonar, face recognition, recommender systems, cloud removal in satelliteimages, anomaly detection, background subtraction for surveillancevideo, system monitoring, etc. Principal component analysis (PCA) is awidely used tool to obtain a low-rank subspace from high dimensionaldata. For a given rank r, PCA finds the r-dimensional linear subspacethat minimizes the square error loss between a data vector and itsprojection onto that subspace. Although PCA is widely used, it is notrobust against outliers. With just one grossly corrupted entry, thelow-rank subspace estimated from classic PCA can be arbitrarily far fromthe true subspace. This shortcoming reduces the value of PCA becauseoutliers are often observed in the modern world's massive data. Forexample, data collected through sensors, cameras, websites, etc. areoften very noisy and contain error entries or outliers.

Robust PCA (RPCA) methods have been investigated to provide goodpractical performance with strong theoretical performance guarantees,but are typically batch algorithms that require loading of allobservations into memory before processing. This makes them inefficientand impractical for use in processing very large datasets commonlyreferred to as “big data” due to the amount of memory required and theslow processing speeds. The robust PCA methods further fail to addressthe problem of slowly or abruptly changing subspace.

SUMMARY

In an example embodiment, a non-transitory computer-readable medium isprovided having stored thereon computer-readable instructions that, whenexecuted by a computing device, cause the computing device to update anestimate of one or more principal components for a next observationvector. An initial observation matrix is defined with a first pluralityof observation vectors. A number of the first plurality of observationvectors is a predefined window length. Each observation vector of thefirst plurality of observation vectors includes a plurality of values.Each value of the plurality of values is associated with a variable todefine a plurality of variables. Each variable of the plurality ofvariables describes a characteristic of a physical object. A principalcomponents decomposition is computed using the defined initialobservation matrix. The principal components decomposition includes asparse noise vectors, a first singular value decomposition vector U, anda second singular value decomposition vector ν for each observationvector of the first plurality of observation vectors. A rank r isdetermined based on the computed principal components decomposition. Anext principal components decomposition is computed for a nextobservation vector using the determined rank r. The computed nextprincipal components decomposition is output for the next observationvector. The output next principal components decomposition is monitoredto determine a status of the physical object.

In another example embodiment, a computing device is provided. Thecomputing device includes, but is not limited to, a processor and anon-transitory computer-readable medium operably coupled to theprocessor. The computer-readable medium has instructions stored thereonthat, when executed by the computing device, cause the computing deviceto update an estimate of one or more principal components for a nextobservation vector.

In yet another example embodiment, a method of updating an estimate ofone or more principal components for a next observation vector isprovided.

Other principal features of the disclosed subject matter will becomeapparent to those skilled in the art upon review of the followingdrawings, the detailed description, and the appended claims.

BRIEF DESCRIPTION OF THE DRAWINGS

Illustrative embodiments of the disclosed subject matter will hereafterbe described referring to the accompanying drawings, wherein likenumerals denote like elements.

FIG. 1 depicts a block diagram of a monitoring device in accordance withan illustrative embodiment.

FIG. 2 depicts a flow diagram illustrating example operations performedby the monitoring device of FIG. 1 in accordance with an illustrativeembodiment.

FIG. 3 depicts a flow diagram illustrating additional example operationsperformed by the monitoring device of FIG. 1 in accordance with anillustrative embodiment.

FIG. 4 depicts a flow diagram illustrating additional example operationsperformed by the monitoring device of FIG. 1 in accordance with anillustrative embodiment.

FIGS. 5A, 5B, and 5C depict a flow diagram illustrating additionalexample operations performed by the monitoring device of FIG. 1 inaccordance with an illustrative embodiment.

FIG. 6 graphically illustrates a plurality of observation vectors inaccordance with an illustrative embodiment.

FIG. 7 provides a comparison between a relative error computed for alow-rank matrix using a known robust principal component analysis (RPCA)algorithm, using the operations of FIGS. 2-4, and using the operationsof FIGS. 3-5 with a slowly changing subspace in accordance with anillustrative embodiment.

FIG. 8 provides a comparison between a relative error computed for asparse noise matrix using the known RPCA algorithm, using the operationsof FIGS. 2-4, and using the operations of FIGS. 3-5 with the slowlychanging subspace in accordance with an illustrative embodiment.

FIG. 9 provides a comparison between a proportion of correctlyidentified elements in the sparse noise matrix using the known RPCAalgorithm, using the operations of FIGS. 2-4, and using the operationsof FIGS. 3-5 with the slowly changing subspace in accordance with anillustrative embodiment.

FIG. 10 provides a comparison as a function of time between the relativeerror computed for the low-rank matrix using the known RPCA algorithm,using the operations of FIGS. 2-4, and using the operations of FIGS. 3-5with a rank of 10, a sparsity probability of 0.01, and the slowlychanging subspace in accordance with an illustrative embodiment.

FIG. 11 provides a comparison as a function of time between the relativeerror computed for the low-rank matrix using the known RPCA algorithm,using the operations of FIGS. 2-4, and using the operations of FIGS. 3-5with a rank of 50, a sparsity probability of 0.01, and the slowlychanging subspace in accordance with an illustrative embodiment.

FIG. 12 provides a comparison between a relative error computed for thelow-rank matrix using the known RPCA algorithm, using the operations ofFIGS. 2-4, and using the operations of FIGS. 3-5 with two abrupt changesin the subspace in accordance with an illustrative embodiment.

FIG. 13 provides a comparison between a relative error computed for thesparse noise matrix using the known RPCA algorithm, using the operationsof FIGS. 2-4, and using the operations of FIGS. 3-5 with the two abruptchanges in the subspace in accordance with an illustrative embodiment.

FIG. 14 provides a comparison between the proportion of correctlyidentified elements in the sparse noise matrix using the known RPCAalgorithm, using the operations of FIGS. 2-4, and using the operationsof FIGS. 3-5 with the two abrupt changes in the subspace in accordancewith an illustrative embodiment.

FIG. 15 provides a comparison as a function of time between the relativeerror computed for the low-rank matrix using the known RPCA algorithm,using the operations of FIGS. 2-4, and using the operations of FIGS. 3-5with a rank of 10, a sparsity probability of 0.01, and the two abruptchanges in the subspace in accordance with an illustrative embodiment.

FIG. 16 provides a comparison as a function of time between the relativeerror computed for the low-rank matrix using the known RPCA algorithm,using the operations of FIGS. 2-4, and using the operations of FIGS. 3-5with a rank of 50, a sparsity probability of 0.01, and the two abruptchanges in the subspace in accordance with an illustrative embodiment.

FIG. 17 provides a comparison as a function of rank between a relativedifference between a detected change point time and a true change pointtime using the operations of FIGS. 3-5 with the sparsity probability of0.01, and the two abrupt changes in the subspace in accordance with anillustrative embodiment.

FIG. 18A shows an image capture from an airport video at a first time inaccordance with an illustrative embodiment.

FIG. 18B shows the low-rank matrix computed from the image capture ofFIG. 18A using the operations of FIGS. 3-5 in accordance with anillustrative embodiment.

FIG. 18C shows the low-rank matrix computed from the image capture ofFIG. 18A using the known RPCA algorithm in accordance with anillustrative embodiment.

FIG. 18D shows the sparse noise matrix computed from the image captureof FIG. 18A using the operations of FIGS. 3-5 in accordance with anillustrative embodiment.

FIG. 18E shows the sparse noise matrix computed from the image captureof FIG. 18A using the known RPCA algorithm in accordance with anillustrative embodiment.

FIG. 19A shows an image capture from the airport video at a second timein accordance with an illustrative embodiment.

FIG. 19B shows the low-rank matrix computed from the image capture ofFIG. 19A using the operations of FIGS. 3-5 in accordance with anillustrative embodiment.

FIG. 19C shows the low-rank matrix computed from the image capture ofFIG. 19A using the known RPCA algorithm in accordance with anillustrative embodiment.

FIG. 19D shows the sparse noise matrix computed from the image captureof FIG. 19A using the operations of FIGS. 3-5 in accordance with anillustrative embodiment.

FIG. 19E shows the sparse noise matrix computed from the image captureof FIG. 19A using the known RPCA algorithm in accordance with anillustrative embodiment.

FIG. 20 shows a comparison of principal components computed using movingwindows RPCA in accordance with an illustrative embodiment.

DETAILED DESCRIPTION

Referring to FIG. 1, a block diagram of a monitoring device 100 is shownin accordance with an illustrative embodiment. Monitoring device 100 mayinclude an input interface 102, an output interface 104, a communicationinterface 106, a non-transitory computer-readable medium 108, aprocessor 110, a decomposition application 122, a dataset 124, a datasetdecomposition description 126, and a dataset change points description128. Fewer, different, and/or additional components may be incorporatedinto monitoring device 100.

Input interface 102 provides an interface for receiving information fromthe user or another device for entry into monitoring device 100 asunderstood by those skilled in the art. Input interface 102 mayinterface with various input technologies including, but not limited to,a keyboard 112, a microphone 113, a mouse 114, a display 116, a trackball, a keypad, one or more buttons, etc. to allow the user to enterinformation into monitoring device 100 or to make selections presentedin a user interface displayed on display 116. The same interface maysupport both input interface 102 and output interface 104. For example,display 116 comprising a touch screen provides a mechanism for userinput and for presentation of output to the user. Monitoring device 100may have one or more input interfaces that use the same or a differentinput interface technology. The input interface technology further maybe accessible by monitoring device 100 through communication interface106.

Output interface 104 provides an interface for outputting informationfor review by a user of monitoring device 100 and/or for use by anotherapplication or device. For example, output interface 104 may interfacewith various output technologies including, but not limited to, display116, a speaker 118, a printer 120, etc. Monitoring device 100 may haveone or more output interfaces that use the same or a different outputinterface technology. The output interface technology further may beaccessible by monitoring device 100 through communication interface 106.

Communication interface 106 provides an interface for receiving andtransmitting data between devices using various protocols, transmissiontechnologies, and media as understood by those skilled in the art.Communication interface 106 may support communication using varioustransmission media that may be wired and/or wireless. Monitoring device100 may have one or more communication interfaces that use the same or adifferent communication interface technology. For example, monitoringdevice 100 may support communication using an Ethernet port, a Bluetoothantenna, a telephone jack, a USB port, etc. Data and messages may betransferred between monitoring device 100 and another computing deviceof a distributed computing system 130 using communication interface 106.

Computer-readable medium 108 is a non-transitory electronic holdingplace or storage for information so the information can be accessed byprocessor 110 as understood by those skilled in the art.Computer-readable medium 108 can include, but is not limited to, anytype of random access memory (RAM), any type of read only memory (ROM),any type of flash memory, etc. such as magnetic storage devices (e.g.,hard disk, floppy disk, magnetic strips, . . . ), optical disks (e.g.,compact disc (CD), digital versatile disc (DVD), . . . ), smart cards,flash memory devices, etc. Monitoring device 100 may have one or morecomputer-readable media that use the same or a different memory mediatechnology. For example, computer-readable medium 108 may includedifferent types of computer-readable media that may be organizedhierarchically to provide efficient access to the data stored therein asunderstood by a person of skill in the art. As an example, a cache maybe implemented in a smaller, faster memory that stores copies of datafrom the most frequently/recently accessed main memory locations toreduce an access latency. Monitoring device 100 also may have one ormore drives that support the loading of a memory media such as a CD,DVD, an external hard drive, etc. One or more external hard drivesfurther may be connected to monitoring device 100 using communicationinterface 106.

Processor 110 executes instructions as understood by those skilled inthe art. The instructions may be carried out by a special purposecomputer, logic circuits, or hardware circuits. Processor 110 may beimplemented in hardware and/or firmware. Processor 110 executes aninstruction, meaning it performs/controls the operations called for bythat instruction. The term “execution” is the process of running anapplication or the carrying out of the operation called for by aninstruction. The instructions may be written using one or moreprogramming language, scripting language, assembly language, etc.Processor 110 operably couples with input interface 102, with outputinterface 104, with communication interface 106, and withcomputer-readable medium 108 to receive, to send, and to processinformation. Processor 110 may retrieve a set of instructions from apermanent memory device and copy the instructions in an executable formto a temporary memory device that is generally some form of RAM.Monitoring device 100 may include a plurality of processors that use thesame or a different processing technology.

Decomposition application 122 performs operations associated withdefining dataset decomposition 126 and/or dataset change pointsdescription 128 from data stored in dataset 124. Dataset decomposition126 may be used to monitor changes in data in dataset 124 to supportvarious data analysis functions as well as provide alert/messagingrelated to the monitored data. Some or all of the operations describedherein may be embodied in decomposition application 122. The operationsmay be implemented using hardware, firmware, software, or anycombination of these methods.

Referring to the example embodiment of FIG. 1, decomposition application122 is implemented in software (comprised of computer-readable and/orcomputer-executable instructions) stored in computer-readable medium 108and accessible by processor 110 for execution of the instructions thatembody the operations of decomposition application 122. Decompositionapplication 122 may be written using one or more programming languages,assembly languages, scripting languages, etc. Decomposition application122 may be integrated with other analytic tools. As an example,decomposition application 122 may be part of an integrated dataanalytics software application and/or software architecture such as thatoffered by SAS Institute Inc. of Cary, N.C., USA. For example,decomposition application 122 may be part of SAS® Enterprise Miner™ orof SAS® Viya™ or of SAS® Visual Data Mining and Machine Learning, all ofwhich are developed and provided by SAS Institute Inc. of Cary, N.C.,USA and may be used to create highly accurate predictive and descriptivemodels based on analysis of vast amounts of data from across anenterprise. Merely for further illustration, decomposition application122 may be implemented using or integrated with one or more SAS softwaretools such as Base SAS, SAS/STATO, SAS® High Performance AnalyticsServer, SAS® LASR™, SAS® In-Database Products, SAS® Scalable PerformanceData Engine, SAS/ORO, SAS/ETSO, SAS® Event Stream Processing, SAS®Inventory Optimization, SAS® Inventory Optimization Workbench, SAS®Visual Analytics, SAS In-Memory Statistics for Hadoop®, SAS® ForecastServer, all of which are developed and provided by SAS Institute Inc. ofCary, N.C., USA. Data mining is applicable in a wide variety ofindustries.

Decomposition application 122 may be integrated with other systemprocessing tools to automatically process data generated as part ofoperation of an enterprise, device, system, facility, etc., to identifyany outliers in the processed data, to monitor changes in the data, totrack movement of an object, to process an image, a video, a sound fileor other data files such as for background subtraction, cloud removal,face recognition, etc., and to provide a warning or alert associatedwith the monitored data using input interface 102, output interface 104,and/or communication interface 106 so that appropriate action can beinitiated in response to changes in the monitored data.

Decomposition application 122 may be implemented as a Web application.For example, decomposition application 122 may be configured to receivehypertext transport protocol (HTTP) responses and to send HTTP requests.The HTTP responses may include web pages such as hypertext markuplanguage (HTML) documents and linked objects generated in response tothe HTTP requests. Each web page may be identified by a uniform resourcelocator (URL) that includes the location or address of the computingdevice that contains the resource to be accessed in addition to thelocation of the resource on that computing device. The type of file orresource depends on the Internet application protocol such as the filetransfer protocol, HTTP, H.323, etc. The file accessed may be a simpletext file, an image file, an audio file, a video file, an executable, acommon gateway interface application, a Java applet, an extensiblemarkup language (XML) file, or any other type of file supported by HTTP.

Dataset 124 may include, for example, a plurality of rows and aplurality of columns. The plurality of rows may be referred to asobservation vectors or records (observations), and the columns may bereferred to as variables. Dataset 124 may be transposed. Dataset 124 mayinclude unsupervised data. The plurality of variables may definemultiple dimensions for each observation vector. An observation vectorm_(i) may include a value for each of the plurality of variablesassociated with the observation i. For example, if dataset 124 includesdata related to operation of a vehicle, the variables may include an oilpressure, a speed, a gear indicator, a gas tank level, a tire pressurefor each tire, an engine temperature, a radiator level, etc. Dataset 124may include data captured as a function of time for one or more physicalobjects.

The data stored in dataset 124 may be generated by and/or captured froma variety of sources including one or more sensors of the same ordifferent type, one or more computing devices, etc. The data stored indataset 124 may be received directly or indirectly from the source andmay or may not be pre-processed in some manner. For example, the datamay be pre-processed using an event stream processor such as the SAS®Event Stream Processing Engine (ESPE), developed and provided by SASInstitute Inc. of Cary, N.C., USA. As a result, the data of dataset 124may be processed as it is streamed through monitoring device 100.

As used herein, the data may include any type of content represented inany computer-readable format such as binary, alphanumeric, numeric,string, markup language, etc. The data may be organized using delimitedfields, such as comma or space separated fields, fixed width fields,using a SAS® dataset, etc. The SAS dataset may be a SAS® file stored ina SAS® library that a SAS® software tool creates and processes. The SASdataset contains data values that may be organized as a table ofobservations (rows) and variables (columns) that can be processed by oneor more SAS software tools.

Dataset 124 may be stored on computer-readable medium 108 or on one ormore computer-readable media of distributed computing system 130 andaccessed by monitoring device 100 using communication interface 106,input interface 102, and/or output interface 104. Data stored in dataset124 may be sensor measurements or signal values captured by a sensor,may be generated or captured in response to occurrence of an event or atransaction, generated by a device such as in response to an interactionby a user with the device, etc. The data stored in dataset 124 mayinclude any type of content represented in any computer-readable formatsuch as binary, alphanumeric, numeric, string, markup language, etc. Thecontent may include textual information, graphical information, imageinformation, audio information, numeric information, etc. that furthermay be encoded using various encoding techniques as understood by aperson of skill in the art. The data stored in dataset 124 may becaptured at different time points periodically, intermittently, when anevent occurs, etc. One or more columns of dataset 124 may include a timeand/or date value.

Dataset 124 may include data captured under normal operating conditionsof the physical object. Dataset 124 may include data captured at a highdata rate such as 200 or more observations per second for one or morephysical objects. For example, data stored in dataset 124 may begenerated as part of the Internet of Things (IoT), where things (e.g.,machines, devices, phones, sensors) can be connected to networks and thedata from these things collected and processed within the things and/orexternal to the things before being stored in dataset 124. For example,the IoT can include sensors in many different devices and types ofdevices, and high value analytics can be applied to identify hiddenrelationships, monitor for system changes, etc. to drive increasedefficiencies. This can apply to both big data analytics and real-timeanalytics. Some of these devices may be referred to as edge devices, andmay involve edge computing circuitry. These devices may provide avariety of stored or generated data, such as network data or dataspecific to the network devices themselves. Again, some data may beprocessed with an ESPE, which may reside in the cloud or in an edgedevice before being stored in dataset 124.

Dataset 124 may be stored using various structures as known to thoseskilled in the art including one or more files of a file system, arelational database, one or more tables of a system of tables, astructured query language database, etc. on monitoring device 100 or ondistributed computing system 130. Monitoring device 100 may coordinateaccess to dataset 124 that is distributed across distributed computingsystem 130 that may include one or more computing devices. For example,dataset 124 may be stored in a cube distributed across a grid ofcomputers as understood by a person of skill in the art. As anotherexample, dataset 124 may be stored in a multi-node Hadoop® cluster. Forinstance, Apache™ Hadoop® is an open-source software framework fordistributed computing supported by the Apache Software Foundation. Asanother example, dataset 124 may be stored in a cloud of computers andaccessed using cloud computing technologies, as understood by a personof skill in the art. The SAS® LASR™ Analytic Server may be used as ananalytic platform to enable multiple users to concurrently access datastored in dataset 124. The SAS® Viya™ open, cloud-ready, in-memoryarchitecture also may be used as an analytic platform to enable multipleusers to concurrently access data stored in dataset 124. Some systemsmay use SAS In-Memory Statistics for Hadoop® to read big data once andanalyze it several times by persisting it in-memory for the entiresession. Some systems may be of other types and configurations.

Bold letters are used herein to denote vectors and matrices. ∥a∥₁ and∥a∥₂ represent the l1-norm and l2-norm of vector a, respectively. For anarbitrary real matrix A, ∥A∥_(F) denotes a Frobenius norm of matrix A,∥A∥₊ denotes a nuclear norm of matrix A (sum of all singular values),and ∥A∥₁ denotes the l1-norm of matrix A, where A is treated as avector.

An objective function for robust principal component analysis (RPCA) maydecompose an observed matrix M (such as data in dataset 124) into alow-rank matrix L and a sparse noise matrix S by solving a principalcomponent pursuit (PCP):

$\begin{matrix}{{{\min\limits_{L,S}{L}_{*}} + {\lambda {S}_{1}}},{{subject}\mspace{14mu} {to}}} & (1) \\{{L + S} = {M.}} & (2)\end{matrix}$

As used herein, T denotes a number of observations in dataset 124, and tis an index to a specific observation vector in dataset 124. m_(t) is anobservation vector of dimension m as in m_(t) is represented by mvariables of the plurality of variables in dataset 124. Data referencedas dataset 124 may be streaming observations m_(t)∈

^(m), t=1, . . . , T, which can be decomposed into two parts asm_(t)=l_(t)+s_(t). The first part l_(t) is a vector from a low-ranksubspace U_(t), and the second part s_(t) is a sparse error with supportsize c_(t), where c_(t) represents a number of nonzero elements ofs_(t), such that c_(t)=Σ_(i=1) ^(m) I_(s) _(t) _([i]≠0). The underlyingsubspace U_(t), may or may not change with time t. When U_(t) does notchange over time, the data in dataset 124 are generated from a stablesubspace. Otherwise, the data in dataset 124 are generated from achanging subspace that may be slowly changing, quickly changing, and/orabruptly changing. s_(t) may satisfy a sparsity assumption, i.e.,c_(t)<<m for all t=1, . . . , T. M_(t)=[m₁, . . . , m_(t)] ∈

^(m×t) denotes a matrix of observation vectors until time t, and M, L, Sdenote M_(t), L_(t), S_(t), respectively. For illustration, FIG. 6 showsa plurality of observation vectors such as a first observation vector m₁600-1, a second observation vector m₂ 600-2, a third observation vectorm₃ 600-3, . . . , and a k^(th) observation vector m_(k) 600-k, whereeach observation vector includes m=[ν₁, ν₂, ν₃, . . . , ν_(m)], where ν₁is a value for a first variable of the plurality of variables, ν₂ is avalue for a second variable of the plurality of variables, ν₃ is a valuefor a third variable of the plurality of variables, . . . , ν_(m) is avalue for an m^(th) variable of the plurality of variables.

As understood by a person of skill in the art, various algorithms havebeen developed to solve RPCA-PCP, including accelerated proximalgradient (APG) (see e.g., Z. Lin, A. Ganesh, J. Wright, L. Wu, M. Chen,and Y. Ma, “Fast convex optimization algorithms for exact recovery of acorrupted low-rank matrix,” Computational Advances in Multi-SensorAdaptive Processing (CAMSAP), vol. 61, 2009.) and augmented Lagrangianmultiplier (ALM) (see e.g., Z. Lin, M. Chen, and Y. Ma, “The augmentedLagrange multiplier method for exact recovery of corrupted low-rankmatrices,” arXiv preprint arXiv:1009.5055, 2010.).

Merely for illustration, in implementing RPCA-PCP-ALM, let S_(τ) denotea shrinkage operator S_(τ)=sgn(x)max(|x|−τ, 0), which is extended tomatrices by applying it to each matrix element. Let

_(τ)(X) denote a singular value thresholding operator given by

_(τ)(X)=US_(τ)(Σ)V⁺, where X=UΣV⁺ is any singular value decomposition.Matrices S₀ and Y₀ may be initialized to zero and a loop of equations(1) to (3) below repeated until a convergence criteria (e.g., number ofiterations tolerance value) specified by a user is satisfied:

L ^((k+1))=

_(μ) ⁻¹ (M−S ^((k))+μ⁻¹ Y ^((k)))  (1)

S ^((k+1))=

_(λμ) ⁻¹ (M−L ^((k+1))+μ⁻¹ Y ^((k)))  (2)

Y ^((k+1)) −Y ^((k))+μ(M−L ^((k+1)) +S ^((k+1))).  (3)

_(μ) ⁻¹ is computed by choosing values for constants μ and λ. A singularvalue decomposition (SVD) is computed of X to determine U, Σ, and V⁺.S_(τ) is applied to Σ with τ=μ⁻¹ to shrink the computed singular values.The singular values smaller than μ⁻¹ are set to zero while U and V⁺ arenot modified.

_(λ82) ⁻¹ is computed using equation (2) by applying S_(τ) with τ=λμ⁻¹to shrink a residual of the decomposition. The process is repeated untilconvergence, where an upper indice (e.g., (k), (k+1)) denote aniteration number for the same observation m_(t).

As another illustration, a stochastic optimization RPCA-PCP(STOC-RPCA-PCP) is described by J. Feng, H. Xu, and S. Yan, “Onlinerobust pca via stochastic optimization,” in Advances in NeuralInformation Processing Systems 26 (C. J. C. Burges, L. Bottou, M.Welling, Z. Ghahramani, and K. Q. Weinberger, eds.), pp. 404-412, CurranAssociates, Inc., 2013. STOC-RPCA-PCP starts by minimizing a lossfunction below:

$\begin{matrix}{{{\min\limits_{L,S}{\frac{1}{2}{{M - L - S}}_{F}^{2}}} + {\lambda_{1}{L}_{*}} + {\lambda_{2}{S}_{1}}},} & (4)\end{matrix}$

where λ₁, μ₂ are tuning parameters defined by a user as understood by aperson of skill in the art. An equivalent form of the nuclear norm canbe defined by the following, which states that the nuclear norm for amatrix L whose rank is upper bounded by r has an equivalent form:

$\begin{matrix}{{L}_{*} = {\inf\limits_{{U \in {\mathbb{R}}^{m \times r}},{V \in {\mathbb{R}}^{n \times r}}}{\left\{ {{{\frac{1}{2}{U}_{F}^{2}} + {\frac{1}{2}{V}_{F}^{2}\text{:}\mspace{14mu} L}} = {UV}} \right\}.}}} & (5)\end{matrix}$

Substituting L with UV and plugging equation (5) into equation (4)results in

$\begin{matrix}{{\min\limits_{{U \in {\mathbb{R}}^{m \times r}},{V \in {\mathbb{R}}^{n \times r}}}{\frac{1}{2}{{M - {UV} - S}}_{F}^{2}}} + {\frac{\lambda_{1}}{2}\left( {{U}_{F}^{2} + {V}_{F}^{2}} \right)} + {\lambda_{2}{S}_{1}}} & (6)\end{matrix}$

where U can be seen as a basis for the low-rank subspace and Vrepresents coefficients of the observations with respect to the basis.STOC-RPCA-PCP minimizes an empirical version of the loss function ofequation (6) and processes one observation vector each time instance.Though t is indicated as a time point other types of monotonicallyincreasing values, preferably with a common difference betweensuccessive values, may be used. Given a finite set of observationvectors M_(t)=[m₁, . . . , m_(t)] ∈

^(m×t), the empirical version of the loss function of equation (6) attime point t can be defined by

${{f_{t}(U)} = {{\frac{1}{t}{\sum\limits_{i = 1}^{t}\; {\left( {m_{i},U} \right)}}} + {\frac{\lambda_{1}}{2\; t}{U}_{F}^{2}}}},$

where the loss function for each observation is defined as

${\left( {m_{i},U} \right)}\overset{\Delta}{=}{{\min\limits_{v,s}{\frac{1}{2}{{m_{i} - {Uv} - s}}_{2}^{2}}} + {\frac{\lambda_{1}}{2}{v}_{2}^{2}} + {\lambda_{2}{{s}_{1}.}}}$

Fixing U as U_(t), ν_(t) and s_(t) can be obtained by solving theoptimization problem:

$\left( {v_{t},s_{t}} \right) = {{\underset{v,s}{argmin}\frac{1}{2}{{m_{t} - {Uv} - s}}_{2}^{2}} + {\frac{\lambda_{1}}{2}{v}_{2}^{2}} + {\lambda_{2}{{s}_{1}.}}}$

Assuming that {ν_(i),s_(i)}_(i=1) ^(t), are known, the basis U_(t) canbe updated by minimizing the following function:

${{g_{t}(U)}\overset{\Delta}{=}{{\frac{1}{t}{\sum\limits_{i = 1}^{t}\; \left( {{\frac{1}{2}{{m_{i} - {Uv}_{i}}}_{2}^{2}} + {\frac{\lambda_{1}}{2}{v_{i}}_{2}^{2}} + {\lambda_{2}{s_{i}}_{1}}} \right)}} + {\frac{\lambda_{1}}{2\; t}{U}_{F}^{2}}}},$

which results in the explicit solution

$U_{t} = {\left\lbrack {\sum\limits_{i = 1}^{t}\; {\left( {m_{i} - s_{i}} \right)v_{i}^{T}}} \right\rbrack \left\lbrack {\left( {\sum\limits_{i = 1}^{t}\; {v_{i}v_{i}^{T}}} \right) + {\lambda_{1}I}} \right\rbrack}^{- 1}$

where ν_(i) ^(T) represents the transpose of ν_(i) and I represents anidentity matrix. Referring to FIG. 2, example operations associated withdecomposition application 122 are described. Decomposition application122 may implement an algorithm referenced herein as online, movingwindow, RPCA-PCP (OMWRPCA). Additional, fewer, or different operationsmay be performed depending on the embodiment of decompositionapplication 122. The order of presentation of the operations of FIG. 2is not intended to be limiting. Although some of the operational flowsare presented in sequence, the various operations may be performed invarious repetitions, concurrently (in parallel, for example, usingthreads and/or distributed computing system 130), and/or in other ordersthan those that are illustrated. A user may execute decompositionapplication 122, which causes presentation of a first user interfacewindow, which may include a plurality of menus and selectors such asdrop down menus, buttons, text boxes, hyperlinks, etc. associated withdecomposition application 122 as understood by a person of skill in theart. The plurality of menus and selectors may be accessed in variousorders. An indicator may indicate one or more user selections from auser interface, one or more data entries into a data field of the userinterface, one or more data items read from computer-readable medium 108or otherwise defined with one or more default values, etc. that arereceived as an input by decomposition application 122.

In an operation 200, one or more indicators may be received thatindicate values for input parameters used to define execution ofdecomposition application 122. As an example, the one or more indicatorsmay be received by decomposition application 122 after selection from auser interface window, after entry by a user into a user interfacewindow, read from a command line, read from a memory location ofcomputer-readable medium 108 or of distributed computing system 130,etc. In an alternative embodiment, one or more input parametersdescribed herein may not be user selectable. For example, a value may beused or a function may be implemented automatically or by default.Illustrative input parameters may include, but are not limited to:

-   -   an indicator of dataset 124, such as a location and a name of        dataset 124;    -   an indicator of a plurality of variables of dataset 124 to        define m, such as a list of column names that include numeric        values;    -   an indicator of a value of T, a number of observations to        process that may be an input or determined by reading dataset        124 or may be determined by incrementing a counter for each        streaming observation processed;    -   an indicator of one or more rules associated with selection of        an observation from the plurality of observations of dataset 124        where a column may be specified with a regular expression to        define the rule;    -   an indicator of a value of n_(win), a window size for a number        of observations (in general, the value of n_(win) should be long        enough to capture a temporal content of the signal, but short        enough that the signal is approximately stationary within the        window; thus, the value is related to how fast the underlying        subspace is changing such that, for an underlying subspace that        changes quickly, a small value is used; for illustration a        default value may be 200);    -   an indicator of a value of n_(start), a number of observations        used to initialize the OMWRPCA algorithm where a default value        may be n_(start)=n_(win) (n_(start) may not be equal to n_(win),        for example, when initialization may include more observations        for accuracy);    -   an indicator of a value of rank r of the r-dimensional linear        subspace that is not used if the n_(start) observations are used        initialize the OMWRPCA algorithm (in general, the value of r may        be selected to be large enough to capture “most” of the variance        of the data, where “most” may be defined, for example, as 90%,        or a maximum allowed value based on a dimensionality of the        observations);    -   an indicator of whether to use the n_(start) observations to        initialize the OMWRPCA algorithm and to define the rank r or to        use zero matrices and the indicated rank r (values of rank r        and/or n_(start) may indicate that one or the other is used as        well. For example, rank r equal to zero may indicate that the        n_(start) observations are used to initialize the OMWRPCA        algorithm without a separate input parameter);    -   an indicator of a value of λ, a first tuning parameter defined        above (in general, the value of λ may be set to either of the        values of λ₁ or λ₂ though other values may be used, for example,        as described in E. J. Candés, X. Li, Y. Ma, and J. Wright,        “Robust principal component analysis?”, Journal of the ACM        (JACM), vol. 58, no. 3, p. 11, 2011);    -   an indicator of a value of μ, a second tuning parameter defined        above (in general, the value of μ may be selected based on a        desired level of accuracy. For illustration,

$\left. {\mu = \frac{{mn}_{win}}{\sum\limits_{i = 1}^{m}\; {\sum\limits_{j = 1}^{n_{win}}\; {M_{ij}}}}} \right);$

-   -   an indicator of a value of λ₁, a third tuning parameter defined        above where a default value may be λ₁=1/√{square root over        (max(m,n_(win)))};    -   an indicator of a value of λ₂, a fourth tuning parameter defined        above where a default value may be λ₂=100/√{square root over        (max(m,n_(win)))};    -   an indicator of an RPCA algorithm to apply such as ALM or APG;    -   an indicator of a maximum number of iterations to perform by the        indicated RPCA algorithm;    -   an indicator of a convergence tolerance to determine when        convergence has occurred perform by the indicated RPCA        algorithm;    -   an indicator of a low-rank decomposition algorithm such as        singular value decomposition (SVD) or principal component        analysis (PCA) (in general, PCA requires computing the        covariance matrix first such that SVD may be appropriate when        computing the covariance matrix is expensive, for example, when        m is large;    -   an indicator of parameters to output such as the low-rank data        matrix for each observation L_(t), the low-rank data matrix        L_(T), the sparse data matrix for each observation S_(t), the        sparse data matrix S_(T), an error matrix that contains noise in        dataset 124, the SVD diagonal vector, the SVD left vector, the        SVD right vector, the value of rank r, a matrix of principal        component loadings, a matrix of principal component scores,        etc.;    -   an indicator of where to store and/or to present the indicated        output such as a name and a location of dataset decomposition        description 126 or an indicator that the output is to a display        in a table and/or is to be graphed, etc.;

The input parameters may vary based on the type of data of dataset 124and may rely on the expertise of a domain expert.

In an operation 201, an iteration index t=1 is initialized.

In an operation 202, a determination is made concerning whether to usethe n_(start) observations to initialize the OMWRPCA algorithm. When then_(start) observations are used to initialize the OMWRPCA algorithm,processing continues in an operation 206. When the n_(start)observations are not used to initialize the OMWRPCA algorithm,processing continues in an operation 204.

In operation 204, the rank r is set to the indicated value, a U₀ matrixis initialized to an (m×r) zero matrix, an A₀ matrix is initialized toan (r×r) zero matrix, a B₀ matrix is initialized to an (m×r) zeromatrix, and processing continues in operation 210.

In operation 206, an initial observation matrix M^(s)=[m_(−(n) _(start)⁻¹⁾, . . . , m₀] is selected from dataset 124, where m is an observationvector having the specified index in dataset 124. The user may or maynot have an initial dataset available to select the initial observationmatrix M^(s). The user also may monitor the first few observations todetermine an appropriate “burn-in period”. The user may specify a firstobservation that is the burn-in number of observations after thebeginning of dataset 124.

In an operation 208, the rank r, the U₀ matrix, the A₀ matrix, and theB₀ matrix are computed, for example, using RPCA-PCP-ALM as shownreferring to FIG. 3 in accordance with an illustrative embodiment.Though RPCA-PCP-ALM is used for illustration, other PCA algorithms maybe used to initialize the parameters described in other manners.

Referring to FIG. 3, example operations associated with decompositionapplication 122 performing RPCA-PCP-ALM are described. Additional,fewer, or different operations may be performed depending on theembodiment of decomposition application 122. The order of presentationof the operations of FIG. 3 is not intended to be limiting. Althoughsome of the operational flows are presented in sequence, the variousoperations may be performed in various repetitions, concurrently (inparallel, for example, using threads and/or distributed computing system130), and/or in other orders than those that are illustrated.

In an operation 300, the S₀ matrix and the Y₀ matrix are initialized to(m×n_(start)) zero matrices, and index k=0 is initialized.

_(μ) ⁻¹ is initialized as

_(τ)(X)=US_(τ)(Σ)V⁺, where X=UΣV⁺ is any singular value decomposition ofM^(s) as described above.

_(λμ) ⁻¹ is initialized by applying S_(τ) with τ=λμ⁻¹ to shrink aresidual of the decomposition.

In an operation 301, increment k using k=k+1.

In an operation 302, L^(k)=

_(μ) ⁻¹ (M^(s)−S^((k−1))+μ⁻¹Y^((k−1))) is computed using the value ofthe second tuning parameter μ.

In an operation 304, S^(k)=

_(λμ) ⁻¹ (M^(s)−L^(k)+μ⁻¹Y^((k−1))) is computed.

In an operation 306, Y^(k)=Y^((k−1))+μ(M^(s)−L^(k)+S^(k)) is computed.

In an operation 308, a determination is made concerning whether or not asolution has converged. When the solution has converged, processingcontinues in an operation 310. When the solution has not converged,processing continues in an operation 301 to repeat operations 301 to308. For example, the maximum number of iterations and/or theconvergence tolerance may be used to determine when convergence hasoccurred. For illustration, convergence may be determined when aFrobenius norm of a residual M^(s)−L^(k)+S^(k) is smaller than apredefined tolerance that may be defined by a user or a default valuemay be used.

In operation 310, the rank r, the U₀ matrix, A_(t), B_(t), s₀ and ν₀ arecomputed from the converged solution for L^(k) and S^(k), and processingcontinues in operation 210, where L^(k)=[l_(−(n) _(start) ⁻¹⁾, . . . ,l₀]=[U_(−(n) _(start) ⁻¹⁾ν_(−(n) _(start) ⁻¹⁾, . . . , U₀ν₀] andS^(k)=[s_(−(n) _(start) ⁻¹⁾, . . . , s₀]. From SVD, L^(k)=Û{circumflexover (Σ)}{circumflex over (V)}, where Û∈

^(m×r), {circumflex over (Σ)}∈

^(r×r), and {circumflex over (V)}∈

^(r×n) ^(start) . U₀=Û{circumflex over (Σ)}^(1/2)∈

^(m×r), A₀=Σ_(i=−(n) _(start) ⁻¹⁾ ⁰ν_(i)ν_(i) ^(T)∈

^(r×r), B₀=Σ_(i=−(n) _(start) ⁻¹⁾ ⁰(m_(i)−s_(i))ν_(i) ^(T)∈

^(m×r). The rank r of the r-dimensional linear subspace is a dimensionof Û, {circumflex over (Σ)}, or {circumflex over (V)}. The ν_(i) vectorsare the columns of {circumflex over (V)} obtained from the SVD of L^(k).

Referring again to FIG. 2, in operation 210, a next observation vectorm_(t), such as m₁ on a first iteration, m₂ on a second iteration, etc.is loaded from dataset 124. For example, the next observation vectorm_(t) is loaded by reading from dataset 124 or by receiving a nextstreamed observation vector.

In an operation 212, a decomposition is computed using the loadedobservation vector m_(t) to update the U_(t) matrix, the A_(t) matrix,the B_(t) matrix, s_(t), and ν_(t), for example, using STOC-RPCA-PCP asshown referring to FIG. 4 modified to use the indicated value of n_(win)in accordance with an illustrative embodiment.

Referring to FIG. 4, example operations associated with decompositionapplication 122 performing RPCA-PCP-ALM are described. Additional,fewer, or different operations may be performed depending on theembodiment of decomposition application 122. The order of presentationof the operations of FIG. 4 is not intended to be limiting.

In an operation 400, ν_(t) and s_(t) are computed by solving theoptimization problem

$\left( {v_{t},s_{t}} \right) = {{\underset{v,s}{argmin}\frac{1}{2}{{m_{t} - {U_{t - 1}v} - s}}_{2}^{2}} + {\frac{\lambda_{1}}{2}{v}_{2}^{2}} + {\lambda_{2}{{s}_{1}.}}}$

In an operation 402, A_(t) is computed using A_(t)=A_(t-1)+ν_(t)ν_(t)^(T)−ν_(t-n) _(win) ν_(t-n) _(win) ^(T), where the value ν_(t-n) _(win)ν_(t-n) _(win) ^(T) is removing the effect of the edge of the window.

In an operation 404, B_(t) is computed usingB_(t)=B_(t-1)+(m_(t)−s_(t))ν_(t) ^(T)−(m_(t-n) _(win) −s_(t-n) _(win))ν_(t-n) _(win) ^(T), where the value (m_(t-n) _(win) −s_(t-n) _(win))ν_(t-n) _(win) ^(T) is removing the effect of the edge of the window.

In an operation 406, Ã is computed using Ã=A_(t)+λ₁I, where A_(t)=[a₁, .. . , a_(r)] ∈

^(r×r), λ₁ is the value of the third tuning parameter, I is an identitymatrix, U_(t-1)=[u₁, . . . u_(r)] ∈

^(m×r), B_(t)=[b₁, . . . , b_(r)] ∈

^(r×r), and Ã=[ã₁, . . . , ã_(r)] ∈

^(r×r). A rank index j is initialized as j=1.

In an operation 408, {tilde over (μ)}₁ is computed using

${\overset{\sim}{u}}_{j} = {{\frac{1}{\overset{\sim}{A}\left\lbrack {j,j} \right\rbrack}\left( {b_{j} - {U_{t - 1}{\overset{\sim}{a}}_{j}}} \right)} + {u_{j}.}}$

In an operation 410, μ_(j) is computed using

$u_{j} = {\frac{1}{\max \left( {{{\overset{\sim}{u}}_{j}}_{2},1} \right)}{{\overset{\sim}{u}}_{j}.}}$

In an operation 412, a determination is made concerning whether or noteach rank has been processed, for example, by comparing the rank index jto the rank r. When each rank has been processed, processing continuesin an operation 416. When each rank has not been processed, processingcontinues in an operation 414.

In operation 414, the rank index j is incremented as j=j+1, andprocessing continues in operation 408 to compute the next values.

In operation 416, U_(t) is updated with the computed values of u_(j),j=1, . . . r, where U_(t)

argmin ½ Tr [U^(T) (A_(t)+λ₁I) U^(T)−Tr(U^(T)B_(t)), and processingcontinues in an operation 214 with the updated values of A_(t), B_(t),U_(t), ν_(t), and S_(t).

Referring again to FIG. 2, in operation 214, the low-rank matrix L_(T),where L_(T)={l₁, l₂, . . . l_(T)}, is updated with values ofl_(t)=U_(t)ν_(t) computed in operation 212, and the sparse noise matrixS_(T), where S_(T)={s₁, s₂, . . . s_(T)} is updated with the value ofs_(t) computed in operation 212.

In an operation 216, window decomposition data is output based on theindicated output selections. For example, one or more of l_(t), s_(t),ν_(t), U_(t), A_(t), B_(t), etc. may be output by storing to datasetdecomposition description 126, by presenting on display 116 in a tableor a graph, by printing to printer 120 in a table or a graph, etc.Dataset decomposition description 126 may be stored on computer-readablemedium 108 or on one or more computer-readable media of distributedcomputing system 130 and accessed by monitoring device 100 usingcommunication interface 106, input interface 102, and/or outputinterface 104. Window decomposition data further may be output bysending to another display or printer of distributed computing system130 and/or published to another computing device of distributedcomputing system 130 for further processing of the data in dataset 124and/or of the window decomposition data. For example, monitoring of thewindow decomposition data may be used to determine whether or not asystem is performing as expected.

In an operation 218, a determination is made concerning whether or notdataset 124 includes one or more unprocessed observation vector. Whendataset 124 includes one or more unprocessed observation vector,processing continues in an operation 220. When dataset 124 does notinclude one or more unprocessed observation vector, processing continuesin an operation 222. When dataset 124 is streamed to monitoring device100, processing may wait until a full observation vector is received.When dataset 124 is not being streamed, the last observation vector mayhave been read from dataset 124. As another option, the value of Tindicated by the user may be less than the number of observationsincluded in dataset 124, in which case, processing continues tooperation 222 when the number of iterations of operation 210 specifiedby the indicated value of T have been completed.

In operation 220, the iteration index t is incremented using t=t+1, andprocessing continues in operation 210 to load the next observationvector and repeat operations 210 to 218.

In operation 222, decomposition data is output based on the indicatedoutput selections. For example, one or more of L_(T), S_(T), U_(T), etc.may be output by storing to dataset decomposition description 126, bypresenting on display 116 in a table or a graph, by printing to printer120 in a table or a graph, etc. Dataset decomposition description 126may be stored on computer-readable medium 108 or on one or morecomputer-readable media of distributed computing system 130 and accessedby monitoring device 100 using communication interface 106, inputinterface 102, and/or output interface 104. Decomposition data furthermay be output by sending to another display or printer of distributedcomputing system 130 and/or published to another computing device ofdistributed computing system 130 for further processing of the data indataset 124 and/or of the decomposition data. For example, monitoring ofthe decomposition data may be used to determine whether or not a systemis performing as expected. Decomposition data and/or windowdecomposition data may be selected for output.

Referring to FIGS. 5A, 5B, and 5C, example operations associated withdecomposition application 122 are described. Decomposition application122 may implement an algorithm referenced herein as OMWRPCA withhypothesis testing (OMWRPCAHT). Additional, fewer, or differentoperations may be performed depending on the embodiment of decompositionapplication 122. The order of presentation of the operations of FIGS.5A, 5B, and 5C is not intended to be limiting.

Similar to operation 200, in an operation 500, one or more indicatorsmay be received that indicate values for input parameters used to defineexecution of decomposition application 122. As an example, the one ormore indicators may be received by decomposition application 122 afterselection from a user interface window, after entry by a user into auser interface window, read from a command line, read from a memorylocation of computer-readable medium 108 or of distributed computingsystem 130, etc. In an alternative embodiment, one or more inputparameters described herein may not be user selectable. For example, avalue may be used or a function may be implemented automatically or bydefault. Illustrative input parameters may include, but are not limitedto:

-   -   the indicator of dataset 124;    -   the indicator of a plurality of variables of dataset 124 to        define m;    -   the indicator of the value of T, the number of observations to        process;    -   the indicator of the one or more rules associated with selection        of an observation from the plurality of observations of dataset        124;    -   the indicator of the value of n_(win), the window size for a        number of observations;    -   the indicator of the value of n_(start), the number of        observations used to initialize the OMWRPCA algorithm;    -   the indicator of the value of rank r of the r-dimensional linear        subspace that is not used if the n_(start) observations are used        initialize the OMWRPCA algorithm;    -   the indicator of whether to use the n_(start) observations to        initialize the OMWRPCA algorithm;    -   the indicator of the value of λ, the first tuning parameter;    -   the indicator of the value of μ, the second tuning parameter;    -   the indicator of the value of λ₁, the third tuning parameter;    -   the indicator of the value of λ₂, the fourth tuning parameter;    -   the indicator of the RPCA algorithm to apply such as ALM or APG;    -   the indicator of the maximum number of iterations to perform by        the indicated RPCA algorithm;    -   the indicator of the convergence tolerance to determine when        convergence has occurred perform by the indicated RPCA        algorithm;    -   the indicator of the low-rank decomposition algorithm;    -   an indicator of n_(stable), a first observation window length to        process for a stable operation of OMWRPCA where a default value        may be n_(stable)=n_(win),    -   an indicator of n_(sample), a second observation window length        to process after n_(stable) observations and before a start of        hypothesis testing where a default value may be        n_(sample)=n_(win)/2;

an indicator of α, an abnormal observation threshold to indicate anabnormal observation where a default value may be α=0.01;

-   -   an indicator of N_(check), a buffer window length for hypothesis        testing, that may be less than n_(win)/2 to avoid missing a        change point, but not too small to generate an inordinate number        of false alarms (for illustration, N_(check)=n_(win)/10 may be a        default value);    -   an indicator of n_(prob), a change point threshold that defines        a number of observations within the buffer window length used to        indicate a change point within a current buffer window where a        default value may be n_(prob)=a_(prob)N_(check), where        α_(prob)=0.5 is a default value though the value can be selected        based on typical rates of change and ranges of variation in the        operating regime of dataset 124, where noisier signals may use        longer windows to average out the noise and fast-changing        systems need shorter windows to capture the change points;    -   an indicator of n_(positive), a number of consecutive abnormal        indicators to identify a change point where n_(positive)=3 is a        default value;    -   an indicator of n_(tol), a tolerance parameter that may have a        default value of zero;    -   the indicator of parameters to output such as the low-rank data        matrix for each observation L_(t), the low-rank data matrix        L_(T), the sparse data matrix for each observation S_(t), the        sparse data matrix S_(T), an error matrix that contains noise in        dataset 124, the SVD diagonal vector, the SVD left vector, the        SVD right vector, a matrix of principal component loadings, a        matrix of principal component scores, a change point index, a        change point value, etc.;    -   the indicator of where to store and/or to present the indicated        output such as a name and a location of dataset decomposition        description 126 or an indicator that the output is to a display        in a table and/or is to be graphed, etc.;    -   the indicator of where to store and/or to present the indicated        output such as a name and a location of dataset change points        description 128 or an indicator that the output is to a display        in a table and/or is to be graphed, etc.;

In an operation 502, iteration indices t=0 and t_(start)=0 areinitialized. Additionally, a zero counter vector H_(c) is initialized toa 0^(m+1) zero vector, a flag buffer list B_(f) is initialized as anempty list, and a zero count buffer list B_(c) is initialized as anempty list.

Similar to operation 202, in an operation 504, a determination is madeconcerning whether to use the n_(start) observations to initialize theOMWRPCA algorithm. When the n_(start) observations are used toinitialize the OMWRPCA algorithm, processing continues in an operation510. When the n_(start) observations are not used to initialize theOMWRPCA algorithm, processing continues in an operation 506.

Similar to operation 204, in operation 506, the rank r is set to theindicated value, a U₀ matrix is initialized to an (m×r) zero matrix, anA₀ matrix is initialized to an (r×r) zero matrix, a B₀ matrix isinitialized to an (m×r) zero matrix, and processing continues inoperation 514.

Similar to operation 206, in operation 510, an initial observationmatrix M^(s)=[m_(−(n) _(start) ⁻¹⁾, . . . , m₀] is selected from dataset124, where m is an observation vector having the specified index indataset 124.

Similar to operation 208, in an operation 512, the rank r, the U₀matrix, the A₀ matrix, and the B₀ matrix are computed, for example,using RPCA-PCP-ALM as shown referring to FIG. 3 in accordance with anillustrative embodiment. Though RPCA-PCP-ALM is used for illustration,other PCA algorithms may be used to initialize the parameters describedin other manners.

Similar to operation 220, in operation 514, the iteration index t isincremented using t=t+1.

Similar to operation 210, in an operation 516, a next observation vectorm_(t), such as m₁ on a first iteration, m₂ on a second iteration, etc.is loaded from dataset 124. For example, the next observation vectorm_(t) is loaded by reading from dataset 124 or by receiving a nextstreamed observation vector.

Similar to operation 212, in an operation 518, a decomposition iscomputed using the loaded observation vector m_(t) to update the U_(t)matrix, A_(t), B_(t), s_(t), and ν_(t), for example, using STOC-RPCA-PCPas shown referring to FIG. 4 modified to use the indicated value ofn_(win) in accordance with an illustrative embodiment.

Similar to operation 214, in an operation 520, the low-rank matrixL_(T), where L_(T)={l₁, l₂, . . . l_(T)}, is updated with values ofl_(t)=U_(t)ν_(t) computed in operation 518, and the sparse noise matrixS_(T), where S_(T)={s₁, s₂, . . . s_(T)} is updated with the value ofs_(t) computed in operation 518.

Similar to operation 216, in an operation 522, window decomposition datais output based on the indicated output selections, and processingcontinues in operation 528 shown referring to FIG. 5B.

Similar to operation 218, in an operation 524, a determination is madeconcerning whether or not dataset 124 includes one or more unprocessedobservation vector. When dataset 124 includes one or more unprocessedobservation vector, processing continues in operation 514 to load thenext observation vector and repeat one or more of operations 514 to 570.When dataset 124 does not include one or more unprocessed observationvector, processing continues in an operation 526.

Similar to operation 222, in operation 526, decomposition data is outputbased on the indicated output selections.

Referring to FIG. 5B, in operation 528, a zero counter time series valueĉ_(t) is computed using ĉ_(t)=Σ_(i=1) ^(m)s_(t)[i].

In an operation 530, a determination is made concerning whether or notthe subspace is stable. When the subspace is stable, processingcontinues in operation 532. When the subspace is not stable, processingcontinues in operation 524. Stability may be based on processing anumber of additional observation vectors m_(t) based on the indicatedfirst observation window length n_(stable). For example, whent≥n_(stable)+t_(start) or when t≥n_(stable)+t_(start), the subspace maybe determined to be stable.

In an operation 532, a determination is made concerning whether or not asufficient hypothesis sample size has been processed. When thesufficient hypothesis sample size has been processed, processingcontinues in operation 536. When the sufficient hypothesis sample sizehas not been processed, processing continues in operation 534. Thesufficient hypothesis sample size may be based on processing a number ofadditional observation vectors m_(t) based on the indicated secondobservation window length n_(sample). For example, whent>n_(stable)+n_(sample)+t_(start) or whent≥n_(stable)+n_(sample)+t_(start), the sufficient sample size has beenprocessed. Hypothesis testing is performed on each subsequentobservation vector m_(t).

In operation 534, zero counter vector H_(c) is updated asH_(c)[ĉ_(t)+1]=H_(c)[ĉ_(t)+1]+1, and processing continue in operation524 to process a next observation vector. Zero counter vector H_(c)maintains a count of a number of times each value of ĉ_(t) occurs in acurrent stable period after the last change point.

In operation 536, a probability value p is computed for ĉ_(t) usingp=Σ_(i=ĉ) _(t) _(−n) _(tol) ₊₁ ^(m+1) H_(c)[i]/Σ_(i=0) ^(m+1) H_(c)[i].Probability value p is a probability that a sparse noise vector includesat least as many nonzero elements as the current observation vectorm_(t).

In an operation 538, the computed probability value p is compared to theabnormal observation threshold α.

In an operation 540, a flag value f_(t) is set to zero or one based onthe comparison. The flag value f_(t) is an indicator of whether or notthe current observation vector m_(t) is abnormal in comparison toprevious observation vectors. For example, f_(t)=1 when p<α, or when p≤αto indicate that the current observation vector m_(t) is abnormal, andotherwise, f_(t)=0. Of course, other flag settings and comparisons maybe used as understood by a person of skill in the art. For example,f_(t)=1 may indicate an abnormal observation.

In an operation 542, a determination is made concerning whether or notflag buffer list B_(f) and zero count buffer list B_(c) are full. Whenthe buffer lists are full, processing continues in an operation 544.When the buffer lists are not full, processing continues in an operation548. For example, flag buffer list B_(f) and zero count buffer listB_(c) may be full when their length is N_(check).

In operation 544, a first value in flag buffer list B_(f) and in zerocount buffer list B_(c) is popped from each list to make room foranother item at an end of each buffer list, where c is the value poppedfrom B_(c).

In an operation 546, the H_(c) vector is updated asH_(c)[c+1]=H_(c)[c+1]+1. H_(c) is a zero counter vector that maintains acount of a number of times each value of ĉ_(t) equals the total numberof nonzero components of the estimated sparse vector in the currentstable period after the last change point.

In operation 548, the flag value f_(t) is appended to an end of flagbuffer list B_(f), and ĉ_(t) is appended to an end of zero count bufferlist B_(c).

In an operation 550, a determination is made concerning whether or notflag buffer list B_(f) and zero count buffer list B_(c) are now full.When the buffer lists are now full, processing continues in an operation552 shown referring to FIG. 5C. When the buffer lists are not now full,processing continues in operation 524. For example, flag buffer listB_(f) and zero count buffer list B_(c) may be full when their length isN_(check).

Referring to FIG. 5C, hypothesis testing is performed. In operation 552,a number of abnormal values n_(abnormal) is computed from flag bufferlist B_(f) using n_(abnormal)=Σ_(i=1) ^(N) ^(check) B_(f)[i].

In an operation 554, the computed number of abnormal values n_(abnormal)is compared to the change point threshold n_(prob).

In an operation 556, a determination is made concerning whether or not achange point has occurred based on the comparison. When the change pointhas occurred, processing continues in an operation 558. When the changepoint has not occurred, processing continues in operation 524. Forexample, when n_(abnormal)>n_(prob), or when n_(abnormal)≥n_(prob), achange point has occurred in the most recent N_(check) number ofobservation vectors.

In an operation 558, a buffer starting index t_(bufstt) is identified bylooping through entries in flag buffer list B_(f) from a start of thelist to an end of the list to identify occurrence of a change point,which is detected as a first instance of the number of consecutiveabnormal indicators (i.e., f_(t)=1) n_(positive) in flag buffer listB_(f). For example, if B_(f) [n1], B_(f) [n1+1], B_(f)[n1+n_(positive)−1] are the first occurrence of the n_(positive)consecutive abnormal flag values, n1 is the buffer starting indext_(bufstt).

In an operation 560, a determination is made concerning whether or not abuffer starting index t_(bufstt) was identified. When the bufferstarting index t_(bufstt) was identified, processing continues in anoperation 562. When the buffer starting index t_(bufstt) was notidentified, processing continues in operation 524.

In operation 562, an overall index t_(ov) is determined for theobservation vector associated with the change point based on a currentvalue of t as t_(ov)=t−N_(check)+t_(bufstt). Typically, only the time ofthe change point is important. However, if the change point vector isneeded, it is stored and can be extracted using overall index t_(ov).

In an operation 564, change point data is output based on the indicatedoutput selections. For example, m_(t) _(ov) , may be output by storingto dataset change point description 128, by presenting on display 116 ina table or a graph, by printing to printer 120 in a table or a graph,etc. Dataset change point description 128 may be stored oncomputer-readable medium 108 or on one or more computer-readable mediaof distributed computing system 130 and accessed by monitoring device100 using communication interface 106, input interface 102, and/oroutput interface 104. Change point data m_(t) _(ov) , further may beoutput by sending to another display or printer of distributed computingsystem 130 and/or published to another computing device of distributedcomputing system 130 for further processing of the data in dataset 124and/or of the decomposition data. For example, monitoring of the changepoint data may be used to determine whether or not a system isperforming as expected. Change point data further may include outputthat is associated with m_(t) _(ov) , such as of l_(t) _(ov) , s_(t)_(ov) , ν_(t) _(ov) , U_(t) _(ov) , A_(t) _(ov) , B_(t) _(ov) , etc. Forexample, an alert may be sent when a change point is identified. Thealert may include a reference to t_(ov) or to m_(t) _(ov) or to l_(t)_(ov) , s_(t) _(ov) , ν_(t) _(ov) , U_(t) _(ov) , A_(t) _(ov) , B_(t)_(ov) , etc.

In an operation 566, t_(ov) may be appended to a change point list.

In an operation 568, iteration index t is reset to t=t_(ov) andt_(start)=t_(ov) are initialized

In an operation 570, the H_(c) vector is initialized to the 0^(m+1) zerovector, the flag buffer list B_(f) is initialized as the empty list, andthe zero count buffer list B_(c) is initialized as the empty list, andprocessing continues in operation 504 of FIG. 5A to restart from theidentified change point vector. Processing may continue as windows ofdata are received or read from dataset 124 until processing is stopped.

Referring to FIG. 7, a comparison between a relative error computed forthe low-rank matrix using three different algorithms is shown inaccordance with an illustrative embodiment. The first algorithm is theknown RPCA algorithm described by FIG. 4 without use of windows orinitialization based on the operations of FIG. 3 also indicated asSTOC-RPCA-PCP. The second algorithm uses the operations of FIGS. 2-4with windowing and initialization based on the operations of FIG. 3 alsoindicated as OMWRPCA. The third algorithm uses the operations of FIGS.3-5 with windowing and initialization based on the operations of FIG. 3and hypothesis testing also indicated as OMWRPCAHT.

The comparison was based on data with a slowly changing subspace. Theobservation vectors in the slowly changing subspace were generatedthrough M=L+S, where S is the sparse noise matrix with a fraction of asparsity probability p of non-zero elements. The non-zero elements of Swere randomly chosen and generated from a uniform distribution over theinterval of [−1000, 1000]. The low-rank subspace L was generated as aproduct L=UV, where the sizes of U and V are m×r and m×T, respectively.The elements of both U and V are independent and identically distributed(i.i.d.) samples from the

(0,1) distribution. U was the basis of the constant subspace withdimension r. T=5000 and m=400 were used to generate the observationvectors M. An n_(start) number of observation vectors with size 400×200was also generated. Four combinations of (r, p) were generated: 1) (10,0.01), (10, 0.1), (50, 0.01), and (50, 0.1). For each combination, 50replications were run using the following parameters for each simulationstudy, λ₁=1/√{square root over (400)}, λ₂=100/√{square root over (400)},n_(win)=200, n_(start)=200, n_(stable)=200, n_(sample)=100,N_(check)=20, n_(prob)=0.5 N_(check), α=0.01, n_(positive)=3, and,n_(tol)=0. Subspace U was changed linearly after generating U₀ by addingnew matrices {Ũ_(k)}_(k=1) ^(K) that were generated independently withi.i.d.

(0,1) elements added to the first r₀ columns of U, where Ũ_(k) ∈

^(m×r) ⁰ , k=1, . . . , K and

$K = {\frac{T}{T_{p}}.}$

T_(p)=250 and r₀=5 were used.

${Specifically},{{{for}\mspace{14mu} t} = {{T_{p}*i} + j}},{{{where}\mspace{14mu} i} = 0},\ldots \mspace{14mu},{K - 1},{j = 0},\ldots \mspace{14mu},{T_{p} - 1},{{U_{t}\left\lbrack {:{,{1\text{:}r_{0}}}} \right\rbrack} = {{U_{0}\left\lbrack {:{,{1\text{:}r_{0}}}} \right\rbrack} + {\Sigma_{1 \leq k \leq i}{\overset{\sim}{U}}_{k}} + {\frac{j}{T_{p}}{\overset{\sim}{U}}_{i + 1}}}},{{U_{t}\left\lbrack {:{,{\left( {r_{0} + 1} \right):r}}} \right\rbrack} = {U_{0}\left\lbrack {:{,{\left( {r_{0} + 1} \right):r}}} \right\rbrack}}$

was used to generate the slowly changing subspace observation vectors.

The comparison between the relative error computed for the low-rankmatrix using the three different algorithms and the slowly changingsubspace is shown for each simulation study using a box plot for each.For example, a first box plot 700 shows results for the first algorithmwith (10, 0.01), where (r=10, p=0.01), a second box plot 702 showsresults for the second algorithm with (10, 0.01), and a third box plot704 shows results for the third algorithm with (10, 0.01). A fourth boxplot 710 shows results for the first algorithm with (10, 0.1), a fifthbox plot 712 shows results for the second algorithm with (10, 0.1), anda sixth box plot 714 shows results for the third algorithm with (10,0.1). A seventh box plot 720 shows results for the first algorithm with(50, 0.01), an eighth box plot 722 shows results for the secondalgorithm with (50, 0.01), and a ninth box plot 724 shows results forthe third algorithm with (50, 0.01). A tenth box plot 730 shows resultsfor the first algorithm with (50, 0.1), an eleventh box plot 732 showsresults for the second algorithm with (50, 0.1), and a twelfth box plot734 shows results for the third algorithm with (50, 0.1).

Referring to FIG. 8, a comparison between a relative error computed forthe sparse noise matrix using the three different algorithms and theslowly changing subspace is shown for each simulation study using a boxplot for each in accordance with an illustrative embodiment. Forexample, a first box plot 800 shows results for the first algorithm with(10, 0.01), where (r=10, p=0.01), a second box plot 802 shows resultsfor the second algorithm with (10, 0.01), and a third box plot 804 showsresults for the third algorithm with (10, 0.01). A fourth box plot 810shows results for the first algorithm with (10, 0.1), a fifth box plot812 shows results for the second algorithm with (10, 0.1), and a sixthbox plot 814 shows results for the third algorithm with (10, 0.1). Aseventh box plot 820 shows results for the first algorithm with (50,0.01), an eighth box plot 822 shows results for the second algorithmwith (50, 0.01), and a ninth box plot 824 shows results for the thirdalgorithm with (50, 0.01). A tenth box plot 830 shows results for thefirst algorithm with (50, 0.1), an eleventh box plot 832 shows resultsfor the second algorithm with (50, 0.1), and a twelfth box plot 834shows results for the third algorithm with (50, 0.1).

Referring to FIG. 9, a comparison between a proportion of correctlyidentified elements in the sparse noise matrix using the three differentalgorithms and the slowly changing subspace is shown for each simulationstudy using a box plot for each in accordance with an illustrativeembodiment. For example, a first box plot 900 shows results for thefirst algorithm with (10, 0.01), where (r=10, p=0.01), a second box plot902 shows results for the second algorithm with (10, 0.01), and a thirdbox plot 904 shows results for the third algorithm with (10, 0.01). Afourth box plot 910 shows results for the first algorithm with (10,0.1), a fifth box plot 912 shows results for the second algorithm with(10, 0.1), and a sixth box plot 914 shows results for the thirdalgorithm with (10, 0.1). A seventh box plot 920 shows results for thefirst algorithm with (50, 0.01), an eighth box plot 922 shows resultsfor the second algorithm with (50, 0.01), and a ninth box plot 924 showsresults for the third algorithm with (50, 0.01). A tenth box plot 930shows results for the first algorithm with (50, 0.1), an eleventh boxplot 932 shows results for the second algorithm with (50, 0.1), and atwelfth box plot 934 shows results for the third algorithm with (50,0.1).

Referring to FIG. 10, a comparison as a function of time between therelative error computed for the low-rank matrix using the threedifferent algorithms with a rank of 10, a sparsity probability of 0.01,and the slowly changing subspace is shown in accordance with anillustrative embodiment. A first curve 1000 shows results for the firstalgorithm, and a second curve 1002 shows results for the secondalgorithm and the third algorithm. As expected, in a slowly changingsubspace, the second and third algorithms perform similarly and muchbetter than the first algorithm that is not responsive to the changingsubspace.

Referring to FIG. 11, a comparison as a function of time between therelative error computed for the low-rank matrix using the threedifferent algorithms with a rank of 50, a sparsity probability of 0.01,and the slowly changing subspace is shown in accordance with anillustrative embodiment. A third curve 1100 shows results for the firstalgorithm, and a fourth curve 1102 shows results for the secondalgorithm and the third algorithm.

Referring to FIG. 12, a comparison between a relative error computed forthe low-rank matrix using the three different algorithms is shown inaccordance with an illustrative embodiment. The comparison was based ondata with two abrupt changes in the subspace at t=1000 and t=2000.T=3000 was used. The same parameters used to create the slowly changingsubspace observation vectors were used to generate the with two abrupt,where the underlying subspaces U were changed and generated completelyindependently for each time window [0,1000], [1000,2000], and[2000,3000].

For example, a first box plot 1200 shows results for the first algorithmwith (10, 0.01), a second box plot 1202 shows results for the secondalgorithm with (10, 0.01), and a third box plot 1204 shows results forthe third algorithm with (10, 0.01). A fourth box plot 1210 showsresults for the first algorithm with (10, 0.1), a fifth box plot 1212shows results for the second algorithm with (10, 0.1), and a sixth boxplot 1214 shows results for the third algorithm with (10, 0.1). Aseventh box plot 1220 shows results for the first algorithm with (50,0.01), an eighth box plot 1222 shows results for the second algorithmwith (50, 0.01), and a ninth box plot 1224 shows results for the thirdalgorithm with (50, 0.01). A tenth box plot 1230 shows results for thefirst algorithm with (50, 0.1), an eleventh box plot 1232 shows resultsfor the second algorithm with (50, 0.1), and a twelfth box plot 1234shows results for the third algorithm with (50, 0.1).

Referring to FIG. 13, a comparison between a relative error computed forthe sparse noise matrix using the three different algorithms and withtwo abrupt changes in the subspace is shown in accordance with anillustrative embodiment. For example, a first box plot 1300 showsresults for the first algorithm with (10, 0.01), where (r=10, p=0.01), asecond box plot 1302 shows results for the second algorithm with (10,0.01), and a third box plot 1304 shows results for the third algorithmwith (10, 0.01). A fourth box plot 1310 shows results for the firstalgorithm with (10, 0.1), a fifth box plot 1312 shows results for thesecond algorithm with (10, 0.1), and a sixth box plot 1314 shows resultsfor the third algorithm with (10, 0.1). A seventh box plot 1320 showsresults for the first algorithm with (50, 0.01), an eighth box plot 1322shows results for the second algorithm with (50, 0.01), and a ninth boxplot 1324 shows results for the third algorithm with (50, 0.01). A tenthbox plot 1330 shows results for the first algorithm with (50, 0.1), aneleventh box plot 1332 shows results for the second algorithm with (50,0.1), and a twelfth box plot 1334 shows results for the third algorithmwith (50, 0.1).

Referring to FIG. 14, a comparison between a proportion of correctlyidentified elements in the sparse noise matrix using the three differentalgorithms and with two abrupt changes in the subspace is shown inaccordance with an illustrative embodiment. For example, a first boxplot 1400 shows results for the first algorithm with (10, 0.01), where(r=10, p=0.01), a second box plot 1402 shows results for the secondalgorithm with (10, 0.01), and a third box plot 1404 shows results forthe third algorithm with (10, 0.01). A fourth box plot 1410 showsresults for the first algorithm with (10, 0.1), a fifth box plot 1412shows results for the second algorithm with (10, 0.1), and a sixth boxplot 1414 shows results for the third algorithm with (10, 0.1). Aseventh box plot 1420 shows results for the first algorithm with (50,0.01), an eighth box plot 1422 shows results for the second algorithmwith (50, 0.01), and a ninth box plot 1424 shows results for the thirdalgorithm with (50, 0.01). A tenth box plot 1430 shows results for thefirst algorithm with (50, 0.1), an eleventh box plot 1432 shows resultsfor the second algorithm with (50, 0.1), and a twelfth box plot 1434shows results for the third algorithm with (50, 0.1).

Referring to FIG. 15, a comparison as a function of time between therelative error computed for the low-rank matrix using the threedifferent algorithms with a rank of 10, a sparsity probability of 0.01,and with two abrupt changes in the subspace is shown in accordance withan illustrative embodiment. A first curve 1500 shows results for thefirst algorithm, a second curve 1502 shows results for the secondalgorithm, and a third curve 1504 shows results for the third algorithm.

Referring to FIG. 16, a comparison as a function of time between therelative error computed for the low-rank matrix using the threedifferent algorithms with a rank of 50, a sparsity probability of 0.01,and with two abrupt changes in the subspace is shown in accordance withan illustrative embodiment. A first curve 1600 shows results for thefirst algorithm, a second curve 1602 shows results for the secondalgorithm, and a third curve 1604 shows results for the third algorithm.

Referring to FIGS. 12-16, as expected, in an abruptly changing subspace,the second algorithm performed much better than the first algorithm thatis not responsive to the changing subspace. The third algorithmperformed much better than the second algorithm that is not responsiveto the abrupt changes.

Referring to FIG. 17, a comparison as a function of rank between arelative difference between a detected change point time and a truechange point time using the third algorithm with the sparsityprobability of 0.01, and the two abrupt changes in the subspace is shownin accordance with an illustrative embodiment. The results show that thethird algorithm can accurately determine the change point time.

Video is a good candidate for low-rank subspace tracking due to acorrelation between frames. In surveillance video, a background isgenerally stable and may change very slowly due to varying illumination.To demonstrate that the third algorithm can effectively track slowlychanging subspace with change points, a “virtual camera” was panned fromleft to right and right to left through an airport video. The “virtualcamera” moved at a speed of 1 pixel per 10 frames. The original framewas size 176×144 for the airport video data. The virtual camera had thesame height and half the width. Each frame was stacked to a column andfed as dataset 124 to the third algorithm. To make the backgroundsubtraction task even more difficult, one change point was added to thevideo where the “virtual camera” jumped instantly from the mostright-hand side to the most left-hand side.

Referring to FIG. 18A, an image capture from an airport video at a firsttime is shown in accordance with an illustrative embodiment. Referringto FIG. 18B, the low-rank matrix computed from the image capture of FIG.18A using the third algorithm is shown in accordance with anillustrative embodiment. Referring to FIG. 18C, the low-rank matrixcomputed from the image capture of FIG. 18A using the first algorithm isshown in accordance with an illustrative embodiment. Referring to FIG.18D, the sparse noise matrix computed from the image capture of FIG. 18Ausing the third algorithm is shown in accordance with an illustrativeembodiment. Referring to FIG. 18E, the sparse noise matrix computed fromthe image capture of FIG. 18A using the first algorithm is shown inaccordance with an illustrative embodiment. The results generated usingthe third algorithm are less blurry with a less populated sparse noisematrix.

Referring to FIG. 19A, an image capture from the airport video at asecond time is shown in accordance with an illustrative embodiment.Referring to FIG. 19B, the low-rank matrix computed from the imagecapture of FIG. 19A using the third algorithm is shown in accordancewith an illustrative embodiment. Referring to FIG. 19C, the low-rankmatrix computed from the image capture of FIG. 19A using the firstalgorithm is shown in accordance with an illustrative embodiment.Referring to FIG. 19D, the sparse noise matrix computed from the imagecapture of FIG. 19A using the third algorithm is shown in accordancewith an illustrative embodiment. Referring to FIG. 19E, the sparse noisematrix computed from the image capture of FIG. 19A using the firstalgorithm is shown in accordance with an illustrative embodiment. Theresults generated using the third algorithm are less blurry with a lesspopulated sparse noise matrix.

Referring to FIG. 20, a comparison of principal components computedusing the third algorithm is shown in accordance with an illustrativeembodiment that shows detection of a wind turbine that begins to deviatefrom the other turbines.

There are applications for decomposition application 122 in areas suchas process control and equipment health monitoring, radar tracking,sonar tracking, face recognition, recommender system, cloud removal insatellite images, anomaly detection, background subtraction forsurveillance video, system monitoring, failure detection in mechanicalsystems, intrusion detection in computer networks, human activityrecognition based on sensor data, etc.

For tracking in video: in each image, the sparse noise matrix Srepresents the pixels of moving objects, while the low-rank subspace Lrepresents a static background. Thus, standard tracking algorithms suchas a Kalman filter may be applied to the sparse noise matrix S onlywhile ignoring the low-rank subspace L.

For recommender systems: the low-rank subspace L may be used torepresent preferences of users according to a rating they assign tomovies, songs or other items. Such preferences can vary with time.Therefore, it is important to monitor change points. Additionally, thesparse noise matrix S may correspond to abnormal/undesirable users suchas bots.

Dataset 124 may include sensor readings measuring multiple key health orprocess parameters at a very high frequency. A change point may indicatewhen a process being monitored deviates from a normal operatingcondition such as when a fault, a disturbance, a state change, etc.occurs. Decomposition application 122 is effective in an onlineenvironment in which high-frequency data is generated.

The word “illustrative” is used herein to mean serving as an example,instance, or illustration. Any aspect or design described herein as“illustrative” is not necessarily to be construed as preferred oradvantageous over other aspects or designs. Further, for the purposes ofthis disclosure and unless otherwise specified, “a” or “an” means “oneor more”. Still further, using “and” or “or” in the detailed descriptionis intended to include “and/or” unless specifically indicated otherwise.

The foregoing description of illustrative embodiments of the disclosedsubject matter has been presented for purposes of illustration and ofdescription. It is not intended to be exhaustive or to limit thedisclosed subject matter to the precise form disclosed, andmodifications and variations are possible in light of the aboveteachings or may be acquired from practice of the disclosed subjectmatter. The embodiments were chosen and described in order to explainthe principles of the disclosed subject matter and as practicalapplications of the disclosed subject matter to enable one skilled inthe art to utilize the disclosed subject matter in various embodimentsand with various modifications as suited to the particular usecontemplated.

What is claimed is:
 1. A non-transitory computer-readable medium havingstored thereon computer-readable instructions that when executed by acomputing device cause the computing device to: define an initialobservation matrix with a first plurality of observation vectors,wherein a number of the first plurality of observation vectors is apredefined window length, wherein each observation vector of the firstplurality of observation vectors includes a plurality of values, whereineach value of the plurality of values is associated with a variable todefine a plurality of variables, wherein each variable of the pluralityof variables describes a characteristic of a physical object; compute aprincipal components decomposition using the defined initial observationmatrix, wherein the principal components decomposition includes a sparsenoise vector s, a first singular value decomposition vector U, and asecond singular value decomposition vector ν for each observation vectorof the first plurality of observation vectors; determine a rank r basedon the computed principal components decomposition; compute a nextprincipal components decomposition for a next observation vector usingthe determined rank r; output the computed next principal componentsdecomposition for the next observation vector; and monitor the outputnext principal components decomposition to determine a status of thephysical object.
 2. The non-transitory computer-readable medium of claim1, wherein the described characteristic of the physical object is apixel of an image of the physical object.
 3. The non-transitorycomputer-readable medium of claim 1, wherein the describedcharacteristic of the physical object is a sensor measurement from asensor connected to detect the characteristic of the physical object. 4.The non-transitory computer-readable medium of claim 3, wherein thesensor is selected from the group consisting of a camera, a videocamera, a radar system, and a sonar system.
 5. The non-transitorycomputer-readable medium of claim 1, wherein the predefined windowlength is associated with a number of observation vectors to initializea robust principal components analysis principal component pursuitcomputation.
 6. The non-transitory computer-readable medium of claim 1,wherein the computed principal components decomposition is computedusing a robust principal components analysis principal component pursuitcomputation that computes a singular value decomposition using all ofthe first plurality of observation vectors.
 7. The non-transitorycomputer-readable medium of claim 1, wherein computing the nextprincipal components decomposition comprises removing an effect of anedge of the predefined window length.
 8. The non-transitorycomputer-readable medium of claim 1, wherein computing the nextprincipal components decomposition for the next observation vectorfurther uses the computed principal components decomposition.
 9. Thenon-transitory computer-readable medium of claim 1, wherein computingthe next principal components decomposition comprises solving anoptimization problem${\left( {v_{t},s_{t}} \right) = {{\underset{v,s}{argmin}\frac{1}{2}{{m_{t} - {U_{t - 1}v} - s}}_{2}^{2}} + {\frac{\lambda_{1}}{2}{v}_{2}^{2}} + {\lambda_{2}{s}_{1}}}},$where m_(t) is the next observation vector, U_(t-1) is a first singularvalue decomposition vector of a most recent observation vector of thefirst plurality of observation vectors, λ₁ is a first predefined tuningparameter, λ₂ is a second predefined tuning parameter, and ∥a∥₁, and∥a∥₂ represent the l1-norm and l2-norm of vector a, respectively. 10.The non-transitory computer-readable medium of claim 9, whereincomputing the next principal components decomposition further comprisesupdating an A_(t) matrix using A_(t)=A_(t-1)+ν_(t)ν_(t) ^(T)−ν_(t-n)_(win) ν_(t-n) _(win) ^(T), where A_(t-1) is an A vector of a mostrecent observation vector of the first plurality of observation vectors,n_(win) is the predefined window length, ν_(t-n) _(win) is a secondsingular value decomposition vector of an oldest observation vector ofthe first plurality of observation vectors, and a^(T) represents atranspose of vector a.
 11. The non-transitory computer-readable mediumof claim 10, wherein computing the next principal componentsdecomposition further comprises updating a B_(t) matrix usingB_(t)=B_(t-1)+(m_(t)−s_(t))ν_(t) ^(T)−(m_(t-n) _(win) −s_(t-n) _(win))ν_(t-n) _(win) ^(T), where B_(t-1) is a B vector of a most recentobservation vector of the first plurality of observation vectors,m_(t-n) _(win) is an oldest observation vector of the first plurality ofobservation vectors, and s_(t-n) _(win) is a sparse noise vector of anoldest observation vector of the first plurality of observation vectors.12. The non-transitory computer-readable medium of claim 11, whereincomputing the next principal components decomposition further comprisesupdating U_(t), a first singular value decomposition vector of the nextobservation vector using U_(t)

argmin ½ Tr[U^(T)(A_(t)+λ₁I)U−Tr(U^(T)B_(t)).
 13. The non-transitorycomputer-readable medium of claim 12, wherein computing the nextprincipal components decomposition further comprises updating L_(t), alow-rank data matrix vector for the next observation vector usingL_(t)=U_(t)ν_(t).
 14. The non-transitory computer-readable medium ofclaim 1, wherein computing the next principal components decompositionis repeated for each of a second plurality of observation vectors as thenext observation vector.
 15. A computing device comprising: a processor;and a non-transitory computer-readable medium operably coupled to theprocessor, the computer-readable medium having computer-readableinstructions stored thereon that, when executed by the processor, causethe computing device to define an initial observation matrix with afirst plurality of observation vectors, wherein a number of the firstplurality of observation vectors is a predefined window length, whereineach observation vector of the first plurality of observation vectorsincludes a plurality of values, wherein each value of the plurality ofvalues is associated with a variable to define a plurality of variables,wherein each variable of the plurality of variables describes acharacteristic of a physical object; compute a principal componentsdecomposition using the defined initial observation matrix, wherein theprincipal components decomposition includes a sparse noise vector s, afirst singular value decomposition vector U, and a second singular valuedecomposition vector ν for each observation vector of the firstplurality of observation vectors; determine a rank r based on thecomputed principal components decomposition; compute a next principalcomponents decomposition for a next observation vector using thedetermined rank r; output the computed next principal componentsdecomposition for the next observation vector; and monitor the outputnext principal components decomposition to determine a status of thephysical object.
 16. The computing device of claim 15, wherein thepredefined window length is associated with a number of observationvectors to initialize a robust principal components analysis principalcomponent pursuit computation.
 17. The computing device of claim 15,wherein computing the next principal components decomposition isrepeated for each of a second plurality of observation vectors as thenext observation vector.
 18. A method of updating an estimate of one ormore principal components for a next observation vector, the methodcomprising: defining, by a computing device, an initial observationmatrix with a first plurality of observation vectors, wherein a numberof the first plurality of observation vectors is a predefined windowlength, wherein each observation vector of the first plurality ofobservation vectors includes a plurality of values, wherein each valueof the plurality of values is associated with a variable to define aplurality of variables, wherein each variable of the plurality ofvariables describes a characteristic of a physical object; computing, bythe computing device, a principal components decomposition using thedefined initial observation matrix, wherein the principal componentsdecomposition includes a sparse noise vector s, a first singular valuedecomposition vector U, and a second singular value decomposition vectorν for each observation vector of the first plurality of observationvectors; determining, by the computing device, a rank r based on thecomputed principal components decomposition; computing, by the computingdevice, a next principal components decomposition for a next observationvector using the determined rank r; outputting, by the computing device,the computed next principal components decomposition for the nextobservation vector; and monitoring, by the computing device, the outputnext principal components decomposition to determine a status of thephysical object.
 19. The method of claim 18, wherein the describedcharacteristic of the physical object is a pixel of an image of thephysical object.
 20. The method of claim 18, wherein the describedcharacteristic of the physical object is a sensor measurement from asensor connected to detect the characteristic of the physical object.21. The method of claim 18, wherein the predefined window length isassociated with a number of observation vectors to initialize a robustprincipal components analysis principal component pursuit computation.22. The method of claim 18, wherein the computed principal componentsdecomposition is computed using a robust principal components analysisprincipal component pursuit computation that computes a singular valuedecomposition using all of the first plurality of observation vectors.23. The method of claim 18, wherein computing the next principalcomponents decomposition comprises removing an effect of an edge of thepredefined window length.
 24. The method of claim 18, wherein computingthe next principal components decomposition for the next observationvector further uses the computed principal components decomposition. 25.The method of claim 18, wherein computing the next principal componentsdecomposition comprises solving an optimization problem${\left( {v_{t},s_{t}} \right) = {{\underset{v,s}{argmin}\frac{1}{2}{{m_{t} - {U_{t - 1}v} - s}}_{2}^{2}} + {\frac{\lambda_{1}}{2}{v}_{2}^{2}} + {\lambda_{2}{s}_{1}}}},$where m_(t) is the next observation vector, U_(t-1) is a first singularvalue decomposition vector of a most recent observation vector of thefirst plurality of observation vectors, λ₁ is a first predefined tuningparameter, λ₂ is a second predefined tuning parameter, and ∥a∥₁ and ∥a∥₂represent the l1-norm and l2-norm of vector a, respectively.
 26. Themethod of claim 25, wherein computing the next principal componentsdecomposition further comprises updating an A_(t) matrix usingA_(t)=A_(t-1)+ν_(t)ν_(t) ^(T)−ν_(t-n) _(win) ν_(t-n) _(win) ^(T), whereA_(t-1) is an A vector of a most recent observation vector of the firstplurality of observation vectors, n_(win) is the predefined windowlength, ν_(t-n) _(win) is a second singular value decomposition vectorof an oldest observation vector of the first plurality of observationvectors, and a^(T) represents a transpose of vector a.
 27. The method ofclaim 26, wherein computing the next principal components decompositionfurther comprises updating a B_(t) matrix usingB_(t)=B_(t-1)+(m_(t)−s_(t))ν_(t) ^(T)−(m_(t-n) _(win) −s_(t-n) _(win))ν_(t-n) _(win) ^(T), where B_(t-1) is a B vector of a most recentobservation vector of the first plurality of observation vectors,m_(t-n) _(win) is an oldest observation vector of the first plurality ofobservation vectors, and s_(t-n) _(win) is a sparse noise vector of anoldest observation vector of the first plurality of observation vectors.28. The method of claim 27, wherein computing the next principalcomponents decomposition further comprises updating U_(t), a firstsingular value decomposition vector of the next observation vector usingU_(t)

argmin ½ Tr[U^(T) (A_(t)+λ₁I)U−Tr(U^(T)B_(t)).
 29. The method of claim28, wherein computing the next principal components decompositionfurther comprises updating L_(t), a low-rank data matrix vector for thenext observation vector using L_(t)=U_(t)ν_(t).
 30. The method of claim18, wherein computing the next principal components decomposition isrepeated for each of a second plurality of observation vectors as thenext observation vector.